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Abstract 

We study the azimuthal modulational instability of vortices with different topological charges, in the focusing two-dimensional 
nonlinear Schrodinger (NLS) equation. The method of studying the stability relies on freezing the radial direction in the Lagrangian 
functional of the NLS in order to form a quasi-one-dimensional azimuthal equation of motion, and then applying a stability analysis 
in Fourier space of the azimuthal modes. We formulate predictions of growth rates of individual modes and find that vortices are 
unstable below a critical azimuthal wave number. Steady state vortex solutions are found by first using a variational approach to 
obtain an asymptotic analytical ansatz, and then using it as an initial condition to a numerical optimization routine. The stability 
analysis predictions are corroborated by direct numerical simulations of the NLS. We briefly show how to extend the method to 
encompass nonlocal nonlinearities that tend to stabilize solutions. 
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1. Introduction 

The nonlinear Schrodinger (NLS) equation has been used 
to describe a very large variety of physical systems since it is 
the lowest order nonlinear (cubic) partial differential equa- 
tion that describes the propagation of modulated waves 
[l] . Two interesting systems described by the NLS that our 
study is relevant to are Bose-Einstein Condensates (BECs) 
|2|.'i| and light propagation in amorphous optical media [3] . 

A BEC is an ultra-cold (on the order of 10~ 8 K) gas of 
10 3 -10 6 atoms which have predominantly condensed into 
the same quantum state, and therefore behaves like one 
large macroscopic atom. Its dynamics can be described 
(through a mean-field approach) by a variant of the NLS 
called the Gross-Pitaevskii (GP) equation that includes an 
external potential trapping the condensed atoms [3] : 

0| *| 2 * + l/ cxt (r)^, (1) 
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where H is the reduced Planck constant, m a is the mass of 
one of the atoms in the condensate, K> x t(r) is the exter- 
nal potential, V 2 is the three-dimensional Laplacian, and 
do is the s-wave scattering length (ao < corresponding 
to the attractive [focusing] case while ao > to the repul- 
sive [defocusing] case). The modulus squared of the wave 
function, I 1 ?! 2 , represents the density of the atoms in the 
condensate. In BECs, a focusing nonlinearity has the phys- 
ical meaning that the particles in the condensate will fea- 
ture attractive interactions. This can cause the BEC to 
collapse into itself, which in turn increases the kinetic en- 
ergy of the particles, and leads to an 'explosive' destruction 
of the BEC dubbed a 'Bosenova' |5|6|7|8j . In the defocus- 
ing case, the particles have repulsive interactions, in which 
case the BEC tries to expand (this is prevented by the ex- 
ternal trap, when the latter is present). Although BECs 
are three-dimensional objects, by increasing the strength 
of the external trap in one transverse direction, one can re- 
shape the BEC into a quasi- two-dimensional disk (or even a 
quasi-one-dimensional cigar-shaped condensate in the case 
of two strong transverse directions) [9] . Each of these sit- 
uations can be described using appropriate forms of the 
two-dimensional (2D) and one-dimensional GP equations 
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On the other hand, amorphous optical media such as sil- 
ica exhibiting a Kerr effect can be modeled using the NLS, 
where the modulus squared of the wave function represents 
the intensity of the light being propagated through the me- 
dia. In such (2 + l)-dimensional NLS is used, where 
the two dimensions of the wave function represent a spa- 
tial cross-section, while the third dimension z (which rep- 
resented time in the case of BECs) represents the direction 
of propagation: 
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where V 2 is the 2D Laplacian, and 0o is the propagation 
constant. The parameters no and n 2 form the index of re- 
fraction in the medium as n — no + 712 1 'I' | 2 , where no is 
the index of refraction in the absence of light, and n 2 is the 
change in the index of refraction due to the intensity of the 
incident light [12] . 

The 2D NLS equation supports vortices [13] . Vortices are 
ring-shaped structures which have a rotational periodic an- 
gular phase associated to them. A key property of the vor- 
tex is its topological charge, denoted as m, which indicates 
how many periods there are in the angular phase around 
the vortex core [3]. For \m\ > 0, the wave function at the 
center of the vortex becomes identically zero, causing the 
ring-like shape. As we will describe below, vortex solutions 
of the NLS in the focusing case are modulationally unsta- 
ble in the azimuthal direction. Our purpose in the present 
manuscript is to formulate and test a method for studying 
the azimuthal modulational instability [14] of vortex solu- 
tions to the NLS. The goal is to predict the growth rates 
of the unstable modes, and predict the critical mode, be- 
low which all modes are unstable. Wherever relevant, we 
will make comparisons of the semi-analytical methods pre- 
sented herein to recent developments in the study of ring 
vortices of the NLS equation, such as Refs. |15I16| . 

The manuscript is organized as follows. In Sec. [2] using 
the Lagrangian representation of the NLS, we formulate a 
quasi-onc-dimcnsional equation of motion for the dynam- 
ics of separable steady-state vortex solutions to the NLS. 
Then, in Sec. [3] we describe the azimuthal modulational 
stability analysis yielding predictions of the growth rates of 
unstable modes, as well as the critical mode, below which 
all modes are unstable. In Sec. 0] a variational approach 
is used to obtain a reliable ansatz for the radial profile of 
steady-state vortices. In Sec. [5] the ansatz is refined into a 
numerically 'exact' radial profile using optimization meth- 
ods and then numerically integrated to extract azimuthal 
growth rates and critical modes that are found to match 
well to our analytical predictions. In Sec. [6] we show an ex- 
tension of the technique to the NLS with nonlocal nonlin- 
earity [17] . Finally, in Sec. [7] we summarize our results and 
give some concluding remarks. 



2. Azimuthal Equation of Motion 

Both physical scenarios described above (BECs and 
amorphous optical media) can be modeled, under appro- 
priate conditions, by the 2D NLS. Let us then use the 
non-dimensionalized NLS 



V 2 * 
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where V 2 ^ is the 2D Laplacian of the wave function W and 
s = +1 (s = —1) denotes the focusing (defocusing) case. 
The action functional of Eq. ([3]) is: 
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Ldt, 
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where the Lagrangian reads 
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L = / CrdrdO, (5) 
Jo Jo 

and its Lagrangian density, in polar coordinates, corre- 
sponds to [18] 

2 



r 



-fl*l 



(0) 



In order to find the azimuthal equation of motion, we 
assume a separable solution with a steady-state "frozen" 
in time radial profile: 

*(r,9,t) = f(r)A(8,t), (7) 

where all of the phase components of the solution are con- 
tained in A, and therefore f(r) £ R. It is worth mentioning 
that vortex solutions to Eq. ([3]) are not necessarily com- 
pletely separable as per Eq. ([7]) and thus this property needs 
to be checked (see Sec. I5.2l for more details). 

When Eq. ([7]) is inserted into Eq. ([5]), since f(r) is 
"frozen" , all radial integrals of Eq. become constants. 
This allows us to transform the 2D Lagrangian into a 
quasi-one dimensional (in 8) Lagrangian which can be 
used to find the equation of motion for A{9, t). We use the 
term 'quasi-one-dimcnsional' because although it becomes 
a one-dimensional problem, the radial direction is not 
ignored, but shows implicitly in the values of the radial 
integral constants. 

First, we insert Eq. ([Jj into the Lagrangian density and 
evaluate the radial integrals of the Lagrangian to obtain 
our quasi-one-dimensional Lagrangian density: 



Ad =^f x {AA\ - A* A t ) + C 2 \A\ 2 + C 3 \A e \ 

+ c 5 a;a + C 6 A e A* - s -c a \a\\ 



(8) 



where 
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We evaluate the variational derivative of the action func- 
tional as shown in Ref. [TH], which in this case takes the 
form: 

SS d dCio d dCio d£m 



5 A* dtd[A* t ] 86 d [At] OA* 



= 0. 



(10) 



Inserting Eq. into Eq. (JTUJ) yields the evolution equation 
for A(9,t): 

idAt = C 2 A-C 3 Aoe + (C 5 - C 6 )A 9 - sd\A\ 2 A. 

Since f(r) is real-valued, and C5 = Cg, the Ag term drops 
out: 

idAt = C 2 A-C 3 A ee -sC 4 \A\ 2 A. (11) 
Applying the rescalings 



-i^rt ) 



and 



C3 



yields the azimuthal NLS 



iA t = -A gg - s-^-\A\ 2 A, 
^3 



(12) 
(13) 

(14) 



that we next study for its modulational instability. 
3. Stability Analysis 

For the stability analysis, we assume a vortex solution of 
Eq. (J3): 

A(0,t) = e l{me+n,t \ (15) 

where m is the topological charge of the vortex, and SI is the 
frequency of rotation of the complex phase. Notice that in 
this context, the vortex waveform becomes an "azimuthal 
plane wave" , and as such its stability analysis becomes the 
standard modulational stability analysis of this plane wave 
(which we briefly review for completeness purposes here) 
|19j . The amplitude of the plane wave does not appear as an 
explicit term because it is absorbed into the radial profile 
f(r) of Eq. Inserting Eq. {TS]) into Eq. ([Tl|). we get the 
following dispersion relation: 

2 , C4 
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Let us now derive equations of motion for a complex 
perturbation. Specifically, we wish to derive the amplitude 
equations for each perturbed Fourier mode. We start by 
perturbing Eq. (|15p with a complex, time-dependent per- 
turbation of the form: 

A(6,t) = (l + u(6,t)+iv(6,t)) e l(me+n,t \ (17) 

where \u\, \v\ <C 1. 

If we rescale time according to the rotating vorticity 
frame as: 

r = t+±-9, 
2m 



this yields 
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s-^- (2uv + u 2 v + v 3 ) 



(18) 
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As in Refs. |19|20j . in order to study modulational insta- 
bility (MI), we seek amplitude equations for the azimuthal 
modes by expanding u and v in a discrete Fourier series: 



M) = ^ £ u(K,t)e 



LKO 



(19) 



K=-c 



where K is the mode number and its respective amplitude 
is given by: 



u(K,t) 



2n 



u{6,t) e tKO d6 



(20) 



v{K,t) = / v{0,t)e 
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Applying these to Eq. (|T8|) yields two coupled nonlinear 
ordinary differential equations describing the dynamics for 
the amplitudes of u and v for each mode. Since we are not 
interested in the long-term dynamics of the system, but 
only in the MI of small perturbations, we drop the nonlinear 
terms and write the resulting linearized system in matrix 
form as: 
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V 



(21) 



The eigenvalues for this linear system are: 



We notice that for the defocusing nonlinearity (s = —1) 
dark vortices are supported by a non-zero background, and 
thus the d integrals do not converge, and therefore the 
method employed above would need to be adjusted by ap- 
propriately subtracting the background field in the La- 
grangian integrals. Nonetheless, it is worth mentioning that 
higher (m > 1) charge dark vortices are unstable since 
they tend to split into unitary charge vortices as shown in 
Ref. [2T] (see also Refs. |22I23I24I25| . and references therein, 
for recent work on this topic). However, this instability is 
not of the modulational type and thus we do not study it 
here, and therefore we concentrate on the focusing case of 
s = +1 ('bright' vortices) below. It is also interesting to 
note that the presence of a confining potential might sta- 
bilize higher order dark vortices in certain parameter win- 
dows |22l26l27ll6j . 
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Returning to the case of interest, namely the focusing 
case (s = +1), there is a bifurcation at a critical value ofK: 



-Kcrit = 



(22) 



where X < Jf cr it indicates a modulational instability. An 
example of such MI is shown in Fig.[T] 




Fig. 1. A typical numerical simulation of a vortex solution to the NLS 
showing MI. The vortex shown is of charge m = 2 perturbed with 
mode K = 5 starting with a perturbation amplitude of e = 0.001. 
(a) t = 0, (b) t = 8, (c) t = 10, and (d) t = 12. 

To predict the actual exponential growth rates for the 
perturbation of each mode from the eigenvalues, the time 
rescaling of Eq. Q13p needs to be taken into account, in 
which case the growth rates (in terms of i^ C rit) are: 

1 'k 2 (Kl, t -K*). (23) 



A± = ± 



4. Variational Approach 

Explicit solutions for two dimensional steady-state vor- 
tices of the NLS are not available. Therefore, in order to find 
a tractable, approximate, solution, we use a variational ap- 
proach (VA) to get a reasonable ansatz, and then use that 
ansatz as an initial condition to a nonlinear equation opti- 
mization routine which finds the numerically 'exact' steady- 
state profile, f(r). The VA-inferred seed may also be of 
value as an initial guess to other numerical techniques that 
have been previously used to obtain such vortices, including 
shooting methods [15] or Newton-type, fixed point schemes 
|16j . The modified Gauss-Newton scheme presented below 
is intended as an alternative to the former ones. 

To perform the VA, we use the technique described in 
Ref. |18j . We insert a vortex ansatz with variable param- 
eters into the Lagrangian of the NLS, and use the Euler- 
Lagrange equations to find the 'best' values for the param- 
eters. We start with a general, separable, steady-state so- 
lution: 

tf(r,0,t) = /(r)e i(m9+nt) , (24) 
where f(r) is the steady-state radial profile which we want 
to find. Inserting this solution into the Lagrangian density 
of the NLS yields: 

2vr ( Q. Ci + m 2 C z + C 2 -\ C 4 ) , (25) 



where we have now explicitly set s = +1 and the C- 
constants are the same as in Eq. ©. 

We use a one-dimensional soliton sech ansatz similar to 
that used in Ref. [28l: 



f(r) = VBsech ( y ^(i — r c 



(26) 



with parameters B and r c corresponding, respectively, to 
the amplitude and location of the ring induced by the vor- 
tex. Now assuming r c to be large, we can approximate the 
C-constants as follows: 
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C 2 = 



/ 5 sech 2 (F)rdr = 2 In (1 + E) 
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r c V8B, (27) 
(28) 
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where E 
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" and F = Jf (r - r c ). The C 3 integral 



does not converge due to the singularity at r = 0. However, 
since we assume r c to be large, the r in the integrand can 
be viewed as a constant (which we choose to be the center 
of the sech, i.e. r c ) and so we have: 

C 3 = / - B sech 2 (F) dr w — / B sech 2 (F) dr 
Jo r r c Jo 

V8B e^ r ' 
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Using these approximations, we obtain: 



(30) 
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The Euler-Lagrange equations: 

dB ' 
lead us to the solution: 

B = 30, 
and so our VA ansatz is: 



2m 2 

IT 



fir) 



'3 O sech 




(31) 



(32) 



Despite the obvious problem with this solution at r = 
(where /(0) should identically be equal to 0), it captures 
the shape and position of the numerically 'exact' solution 
very well (see Fig. Also for higher m, its value at r = 
becomes very close to zero. 

Using the VA ansatz with the asymptotic approximations 
of Eqs. (|27|) - ([30"|) . we can calculate analytical expressions 
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for the exponential growth rates and critical modes of the 
MI: 



K^ t = ±2V2m, 



(33) 



V8m 2 - K 2 
2^ ' 



The advantage of the above formula is that, although ap- 
proximate, they describe in simple terms the MI experi- 
enced by the vortex. Also, it should be noted that this an- 
alytical prediction becomes more accurate for higher order 
vortices as the VA is able to closely match the actual solu- 
tion as depicted in Fig. [5J 

5. Numerical Results 

5.1. Numerical Optimization 

To refine the ansatz profile into a numerically 'exact' 
solution, we implement a nonlinear optimization scheme 
based on a modified Gauss- Newton scheme [29]. First, we 
insert the following separable steady-state solution into 
Eq. ©: 



*(r,0,t) = fir) 



J(m8+Slt) 



(34) 



which produces an ordinary differential equation for f(r) 
which can be discretized as: 

' fi + D(fi)+ff =0, (35) 



where 



n 



D(fi) = 



1 1 

n Ar 



fi+l — fi 



fi fi 

Ar 



(36) 
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We now want a profile, / , which optimizes F towards the 
specific value 0. To do this we iterate the trial profile using: 

fk+i = fk+ a kP k , 
where the step length a&, is found by: 

a Pk) -* a k, 

and where M{f k ) is the merit function defined by: 

n 

M{f) = -Y,mf)) 2 - (37) 

i=i 

The step direction, p k is found using a modified Gauss- 
Newton (GN) formulation: 

Pk = -{JlJ k + Akir 1 JlF(f k ), (38) 

where A^. is called the forcing term, which ensures that the 
step is always defined, even near non-zero roots of M . The 
forcing term is manually set to values which produce desired 
results for our problem (A k = 0.001). Some sample profiles 
for various charges are shown in Fig. [21 where we can see 
the very good agreement between VA and the numerically 
'exact' solution, especially for higher charges. 

The apparent convergence of the VA ansatz with the GN 
refined profile as \m\ increases can be very useful. For very 




Fig. 2. (Color online) Comparison between the VA ansatz 
(dashed/red lines) and the numerically 'exact' solution (solid/blue 
lines) for charges m = 1, ...,7 (curves left to right). We notice that 
the VA captures the numerically 'exact' solution very well as m in- 
creases. 

large \m\, the GN computation using a high enough reso- 
lution to avoid numerical errors can become very compu- 
tationally expensive. Therefore, the analytic stability pre- 
dictions of the VA [see Eq. ([3"3"j) ] can be used for predictions 
without the need to run numerical computations at all. 
Even for low charges, the VA ansatz accurately describes 
the radius and maximum intensity of the vortex. 




Fig. 3. (Color online) Numerical predictions of growth rates of per- 
turbations of azimuthal modes (K) for vortices with Q = 0.25 and 
charges m = 1,...,5 (left to right) using the numerical routine de- 
scribed in Sec.[5]to converge the VA ansatz into a numerically 'exact' 
solution. The predictions are made numerically integrating the con- 
stants of Eq. (0. We see that for each m, after the critical mode, 
the growth rate predictions for each K become indicating that the 
perturbations after the critical mode are stable. 



5.2. Two-dimensional Simulations 

We now compare our predictions for the MI growth rates 
for vortex charges m = 1, 5 using Eq. (|23|) to numerical 
results, see Fig. [3] To verify our predictions we use a polar- 
grid finite-difference scheme where we treat the time deriva- 
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tive separately from the spatial derivatives. For the time 
derivatives, we use the fourth order Runge-Kutta method. 
For the spatial derivatives we use a second-order central 
difference scheme: 



1 # 



2* 



h3 



A9 2 



For our simulations we use f2 = 0.25, Ar = 0.35, AO = 
27r/80, At — 0.0005, with a length of the simulation t max = 
15, and a perturbation amplitude e = 0.00001. 

Using this scheme, along with Dirichlct boundary condi- 
tions, yields the results in Figs. 0] and [B] for m = 2 and m = 
3, respectively. The growth rates are calculated by record- 
ing the maximum and minimum of the modulus squared of 
the crest of the vortex, and computing the average growth 
rate of the perturbation growth. 




Fig. 4. (Color online) Average growth rates from full two-dimen- 
sional simulation of vortices of charge m = 2 perturbed with modes 
K = 1, ...,7 compared to numerical predictions. The predicted 
growth rates are shown in blue circles, while the green squares rep- 
resent the computed growth rates from the full simulation. 
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Fig. 5. (Color online) Same as in Fig. [4] for m = 3. 

Overall, we see that our numerical simulations yield 
growth rates that are close to those predicted, but typi- 
cally slightly higher, with an error on the order of 10% for 
modes far from -Kent- For modes close to -Kcrit, we observe 
higher error. Also, for m — 2 and m = 3, the predicted 
.Kcrit is one mode off. 



Through one-dimensional simulations, as well as numeri- 
cal error analysis, we have accounted for much of this error. 
It is observed that for modes closer to -K C rit, the simulations 
are very sensitive to resolution. By increasing the resolution 
to very high levels, the discrepancy in the one-dimensional 
runs were virtually eliminated. Such high resolutions were 
not used for the 2D simulations because the simulations be- 
come very computationally expensive. Additionally, due to 
the singularities in the C-constant integrals, the numerical 
predictions derived from them also induce slight errors. 

Another source of the discrepancy between our predic- 
tions and the simulations (especially the fact that our criti- 
cal mode prediction is off by one) is that the assumption of 
separability used in Eq. |(7j) is not exact, but rather a good 
approximation. This can be seen by plotting the 2D eigen- 
vectors of the steady-state vortices as seen in Fig. [6] We see 
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Fig. 6. (Color online) Depiction of the modulus squared of numeri- 
cally derived unstable eigenmodes of vortices in the 2D focusing NLS 
of charges m = 1,2,3, and 6 for modes K = 1,...,5 (the vortex of 
charge m = 1 does not have unstable modes past K = 3). It is obvi- 
ous from the panels that the eigenmodes are not completely separa- 
ble into radial and azimuthal parts as assumed in Eq. J7J, but can be 
reasonably approximated by such a separable solution. We also see 
that for higher charges and higher mode numbers, the eigenmodes 
appear to become more separable, and thus the approximation of a 
separable solution becomes more accurate. 

that for low vortex charges, and small mode perturbations, 
the eigenvectors are clearly non-separable into radial and 
azimuthal parts. As one increases the charge and/or the 
mode being perturbed, the eigenvector becomes more sep- 
arable. Since our simulations were done on vortices of low 
charge, some discrepancy due to the assumption of Eq. |(7J) 
is to be expected [3D] . 

Finally, we note in passing that our approach is somewhat 
complementary to the theoretical approach of Ref. [IB] , 
while the combination of both is in some sense tantamount 
to the theoretical analog of what is found numerically in 
Refs. [15|16j . In Ref. [H], the so-called Vakhitov-Kolokolov 
criterion was considered which is implicitly connected to 
the K = perturbation mode and the instability along that 
eigendirection leads to collapse. On the other hand, here 
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we examine the modulational-type instability of higher K 
modes which initiates the unstable dynamics by breakup 
of the azimuthal symmetry (and may, however, eventually 
lead to collapse in conjunction with the K = mode, as 
shown in Fig. [I]). 

6. Nonlocal Nonlinearity 

Here we briefly describe one of the possible extensions 
to the theory, that of incorporating nonlocal interactions. 
Such interactions correspond to various physical systems, 
such as dipole-dipole interactions in a BEC of degenerate 
dipolar atoms 31J, and nonlinear crystals whose nonlinear 
refractive index changes due to the intensity of the light 
present (determined by a transport process such as heat 
conduction) [32] . As we will show, the nonlocality of the 
nonlinearity will induce a stabilizing effect on the mod- 
ulational stability of vortices. Other interesting effects of 
the nonlocality of the nonlinearity include: changing, un- 
der appropriate circumstances, the character of interaction 
of dark solitons from repulsive to attractive [33] ; changing 
the interaction strength between solitons [34]; and stabi- 
lization of dipolc solitons [33J or 2D ring vortices such as 
the ones considered herein [55]. We note that prior work 
has demonstrated that for the case of nonlocal x*- 3 ^ nonlin- 
earity, all three dimensional spatiotcmporal solitons with 
vorticity are unstable [37j . 

For nonlocal interactions, the NLS can be altered to have 
a nonlocal nonlinearity |32] 



V 2 * 



sN \W)y = 



(39) 



where the nonlocal nonlinearity takes the form of a convo- 
lution integral: 

/"2-7T />OC 

N = / V(r' -r,9' -6)\$(r',6',t)\ 2 r'dr'd6', 
Jo Jo 

and where V, the nonlocal response function, is taken to 
be a Gaussian (which appears in relation to the nonlinear 
crystal heat diffusion nonlocality [35] ) : 

f 



V(r' - r, 9' -9) = — - exp 



no' 

where r = (r cos 9, r sin 9) and r' = (r' cos 9', r' s'm9'). 
Formulating the Lagrangian density of Eq. ([39]) yields: 

2 



r 



which, by the same method as in Sec. [5] yields the following 
azimuthal equation of motion: 



iC 1 A t = C 2 A~C 3 A 
where C is defined as: 

/•27T fOO POO 

C(6,t)= / / / 
Jo Jo Jo 

2 tf m l\2 



sC(9,t)A, (40) 
V(r' -r,6' -9) x (41) 



Note that the new, nonlocal, azimuthal NLS ([40]) is the 
same as in the local case [see Eq. ([TT]) ] where the (local) C4 
integral has been replaced by the (nonlocal) convolution C 
integral (|4T]) . 

We use the same stability analysis technique as in Sec. [2] 
with a slight alteration. We notice that if we define: 

/>oo />oo 

R{9'-9)= / V(r' -r,9' -9) f{r) 2 f{r') 2 rr' dr dr' , 
Jo Jo 

then C is now a convolution term as follows: 

C{9,t)= f R{9' -9)\A{9',t)\ 2 d9' = R*\A\ 2 . 



Inserting this into Eq. ([40]) , and using the same rescalings 
as in Eqs. {TSJ) and ([Tg]). yields 

iA t = -Age - ^A(R * \A\ 2 ). 

Now, we perform a stability analysis identical to that of 
Sec. [3J but along with the transforms of u and v, we also 
add: 



2ji 



JKS 



R{K) = / R(9)e l 
Jo 

in which case the convolution term becomes a product 
(AR), and we get: 



x± = ± c[\ 



and, therefore, the critical mode is: 



±\2a 



R(K) 



f(r) 2 f{r') 2 \A(9', t)\ 2 rr' dr dr 1 d& . 



which has the same form as the local nonlinearity case af- 
ter replacing C4 with R(K) . We see that depending on the 
nonlocal response function, K cv i t can be damped, and if 
R(K) < \ then K cl n < 1 and all modes become sta- 
ble. Therefore, we see that one could have a vortex in the 
focusing NLS with a nonlocal nonlinearity which would be 
modulationally stable. In fact, this modulational stability 
(as well as the stability against collapse) of the focusing 
ring vortices of the nonlocal NLS equation has been con- 
firmed in the work of Ref. [31] and is a feature that could 
have practical applications, such as data storage and com- 
munications using light vortices in Kerr optical media |38j . 



7. Conclusions 

We have formulated a methodology for studying az- 
imuthal modulational instability of vortices in the two- 
dimensional nonlinear Schrodingcr (NLS) equation which 
can be extended to incorporate any additional terms in 
the NLS as long as they have a Lagrangian representa- 
tion. (This expandability of the method adds greatly to 
its usefulness and broad relevance). The method relies on 
approximating a vortex solution as being separable into 
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its radial and azimuthal parts, and using the Lagrangian 
functional of the NLS to obtain a quasi-one-dimensional 
equation of motion for the azimuthal direction. A stability 
analysis on modulational perturbations of the equation, 
leads to predictions of growth rates for each perturbed 
mode, and of the critical mode. After obtaining a steady- 
state vortex solution using a variational ansatz along with 
a nonlinear optimization routine, we ran numerical sim- 
ulations of the NLS, perturbing individual modes and 
recording their growth rates to confirm the predictions. 

One key result that should not be overlooked is that of 
the usefulness of the variational ansatz of the vortex pro- 
files that we derived. Since this profile seems to converge 
to the numerically exact solution as the vortex charge be- 
comes large, experimenters can use it to make simple yet ac- 
curate predictions of the vortex radius and intensity given 
experimental parameters. Furthermore, it can be used as 
in Ref . [TS] , in both local and nonlocal settings to yield an 
approximate threshold for collapse dynamics. 

We have also shown theoretical predictions of modula- 
tional instability of vortices which exhibit a nonlocal re- 
sponse by extending the NLS to incorporate a nonlocal non- 
linearity. The results illustrate that nonlocality can damp, 
or completely eliminate, the modulational instability, po- 
tentially leading to the complete stabilization of the non- 
local vortices, as shown numerically, e.g., in Ref. |36j . 
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